#Alexander F. Gazmararian
#afg2@princeton.edu
#January 19, 2024

#Load packages
library(tidyverse)
library(sandwich)
library(lmtest)
library(modelsummary)
library(kableExtra)
library(broom)
library(scales)
library(viridis)
library(here)

#Load custom functions
source(here("code", "fun", "savefig.r"))
source(here("code", "fun", "gofmap.r"))
source(here("code", "fun", "fix_txt.r"))

#Load data
g <- readRDS(here("data", "qualtrics", "natsurvey.rds"))

# Specify empirical model
f.base <- y ~ age + female + black + hispanic + aapi + rural_bin + college + fullemploy + pid3 + ownhome + income5 + buyself_bin + suitable_bin + index_trust + altruism + risk

#Who are the people who don't uptake the treatment
lm(update(f.base, I(worry_num %in% c(3,4)) ~ .), g, subset = (loss == 1)) %>%
  coeftest(., vcovHC(., "HC2"))
#Subset to people who up took the treatment
g.tot <- subset(g, !(worry_num %in% c(3,4) & loss == 1))
#Check share of people who uptook the treatment
nrow(subset(g, loss == 1 & worry_num %in% c(3,4))) / nrow(subset(g, loss == 1))

# Estimate donation
m <- list()
# Donation allocation outcome
m[["Own Donation"]] <- lm(update(f.base, scale(donateloss) ~ loss + .), g.tot)
# Neighbors allocation outcome
m[["Neighbor Donation"]] <- lm(update(f.base, scale(neighbors) ~ loss + .), g.tot)
# Interest in participating
m[["Participate"]] <- lm(update(f.base, interest_num ~ loss + .), g.tot)
# Interest in receiving information
m[["Seek Information"]] <- lm(update(f.base, info_num ~ loss + .), g.tot)
# Share information
m[["Share Information"]] <- lm(update(f.base, share_num ~ loss + .), g.tot)
# Support for program
m[["Support Program"]] <- lm(update(f.base, losssupport_num ~ loss + .), g.tot)

# Create plot data
plot.df <- modelplot(m, vcov = "HC2", draw = FALSE)
plot.df <- subset(plot.df, term == "loss1")
p.tot <-
  plot.df %>%
  mutate(
    lb90 = estimate + std.error * qnorm(.05),
    ub90 = estimate + std.error * qnorm(.95)
  ) %>%
  ggplot(aes(x = estimate, y = model)) +
  geom_vline(xintercept = 0, color = "grey") +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0) +
  geom_errorbar(aes(xmin = lb90, xmax = ub90), width = 0, linewidth = 1.25) +
  geom_point(size = 3, fill = "white", shape = 21) +
  theme_bw(base_size = 13) +
  scale_x_continuous(limits = c(-.1, .35)) +
  scale_y_discrete(limits = rev) +
  labs(y = "Outcome", x = "Complier Average Causal Effect") +
  theme(
    panel.grid = element_blank()
  )
p.tot
savefig(p.tot, "fig_lossaversion", height = 3)

# Heterogeneous treatment effects by party
coefmap <- c(
  "loss1"="Protection Against Loss",
   "loss1:risk"="Protection Against Loss x Risk Preferences",
   "loss1:altruism"="Protection Against Loss x Altruism",
   "loss1:index_trust" = "Protection Against Loss x Trust Index",
   "loss1:buyself_bin"="Protection Against Loss x Solar Interest",
   "loss1:pid3Republican"="Protection Against Loss x Republican",
   "loss1:pid3Neither"="Protection Against Loss x Neither"
  )


# Estimate donation
m.pid <- list()
# Donation allocation outcome
m.pid[["Own"]] <- lm(update(f.base, scale(donateloss) ~ loss * pid3 + .), g.tot)
# Neighbors allocation outcome
m.pid[["Neighbor"]] <- lm(update(f.base, scale(neighbors) ~ loss * pid3 + .), g.tot)
# Interest in participating
m.pid[["Participate"]] <- lm(update(f.base, interest_num ~ loss * pid3 + .), g.tot)
# Support for program
m.pid[["Support"]] <- lm(update(f.base, losssupport_num ~ loss * pid3 + .), g.tot)
# Interest in receiving information
m.pid[["Seek"]] <- lm(update(f.base, info_num ~ loss * pid3 + .), g.tot)
# Share information
m.pid[["Share"]] <- lm(update(f.base, share_num ~ loss * pid3 + .), g.tot)

pid.file <- here("output", "tables", "tab_loss_pid.txt")
modelsummary(
  m.pid,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefmap,
  gof_omit = "R2|IC|Log|F|RMSE|Err",
  add_rows = data.frame(
    c("Covariates"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes")
  ),
  gof_map = gof_map,
  fmt = 2,
  output = "latex"
) %>%
  add_header_above(c(" " = 1, "Donation" = 2, "Program" = 2, "Information" = 2)) %>%
  cat(., file = pid.file)
fix_txt(pid.file)


# Heterogeneous treatment effects by interest in solar
m.interest <- list()
# Donation allocation outcome
m.interest[["Own"]] <- lm(update(f.base, scale(donateloss) ~ loss * buyself_bin + .), g.tot)
# Neighbors allocation outcome
m.interest[["Neighbor"]] <- lm(update(f.base, scale(neighbors) ~ loss * buyself_bin + .), g.tot)
# Interest in participating
m.interest[["Participate"]] <- lm(update(f.base, interest_num ~ loss * buyself_bin + .), g.tot)
# Support for program
m.interest[["Support"]] <- lm(update(f.base, losssupport_num ~ loss * buyself_bin + .), g.tot)
# Interest in receiving information
m.interest[["Seek"]] <- lm(update(f.base, info_num ~ loss * buyself_bin + .), g.tot)
# Share information
m.interest[["Share"]] <- lm(update(f.base, share_num ~ loss * buyself_bin + .), g.tot)

interest.file <- here("output", "tables", "tab_loss_interest.txt")
modelsummary(
  m.interest,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefmap,
  gof_omit = "R2|IC|Log|F|RMSE|Err",
  add_rows = data.frame(
    c("Covariates"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes")
  ),
  gof_map = gof_map,
  fmt = 2,
  output = "latex"
) %>%
  add_header_above(c(" " = 1, "Donation" = 2, "Program" = 2, "Information" = 2)) %>%
  cat(., file = interest.file)
fix_txt(interest.file)


# Heterogeneous treatment effects by risk preferences
m.risk <- list()
# Donation allocation outcome
m.risk[["Own"]] <- lm(update(f.base, scale(donateloss) ~ loss * risk + .), g.tot)
# Neighbors allocation outcome
m.risk[["Neighbor"]] <- lm(update(f.base, scale(neighbors) ~ loss * risk + .), g.tot)
# Interest in participating
m.risk[["Participate"]] <- lm(update(f.base, interest_num ~ loss * risk + .), g.tot)
# Support for program
m.risk[["Support"]] <- lm(update(f.base, losssupport_num ~ loss * risk + .), g.tot)
# Interest in receiving information
m.risk[["Seek"]] <- lm(update(f.base, info_num ~ loss * risk + .), g.tot)
# Share information
m.risk[["Share"]] <- lm(update(f.base, share_num ~ loss * risk + .), g.tot)

risk.file <- here("output", "tables", "tab_loss_risk.txt")
modelsummary(
  m.risk,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefmap,
  gof_omit = "R2|IC|Log|F|RMSE|Err",
  add_rows = data.frame(
    c("Covariates"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes")
  ),
  gof_map = gof_map,
  fmt = 2,
  output = "latex"
) %>%
  add_header_above(c(" " = 1, "Donation" = 2, "Program" = 2, "Information" = 2)) %>%
  cat(., file = risk.file)
fix_txt(risk.file)


# Heterogeneous treatment effects by social preferences
m.altruism <- list()
# Donation allocation outcome
m.altruism[["Own"]] <- lm(update(f.base, scale(donateloss) ~ loss * altruism + .), g.tot)
# Neighbors allocation outcome
m.altruism[["Neighbor"]] <- lm(update(f.base, scale(neighbors) ~ loss * altruism + .), g.tot)
# Interest in participating
m.altruism[["Participate"]] <- lm(update(f.base, interest_num ~ loss * altruism + .), g.tot)
# Support for program
m.altruism[["Support"]] <- lm(update(f.base, losssupport_num ~ loss * altruism + .), g.tot)
# Interest in receiving information
m.altruism[["Seek"]] <- lm(update(f.base, info_num ~ loss * altruism + .), g.tot)
# Share information
m.altruism[["Share"]] <- lm(update(f.base, share_num ~ loss * altruism + .), g.tot)

altruism.file <- here("output", "tables", "tab_loss_altruism.txt")
modelsummary(
  m.altruism,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefmap,
  gof_omit = "R2|IC|Log|F|RMSE|Err",
  add_rows = data.frame(
    c("Covariates"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes")
  ),
  gof_map = gof_map,
  fmt = 2,
  output = "latex"
) %>%
  add_header_above(c(" " = 1, "Donation" = 2, "Program" = 2, "Information" = 2)) %>%
  cat(., file = altruism.file)
fix_txt(altruism.file)


# Heterogeneous treatment effects by risk preferences
m.risk <- list()
# Donation allocation outcome
m.risk[["Own"]] <- lm(update(f.base, scale(donateloss) ~ loss * risk + .), g.tot)
# Neighbors allocation outcome
m.risk[["Neighbor"]] <- lm(update(f.base, scale(neighbors) ~ loss * risk + .), g.tot)
# Interest in participating
m.risk[["Participate"]] <- lm(update(f.base, interest_num ~ loss * risk + .), g.tot)
# Support for program
m.risk[["Support"]] <- lm(update(f.base, losssupport_num ~ loss * risk + .), g.tot)
# Interest in receiving information
m.risk[["Seek"]] <- lm(update(f.base, info_num ~ loss * risk + .), g.tot)
# Share information
m.risk[["Share"]] <- lm(update(f.base, share_num ~ loss * risk + .), g.tot)

risk.file <- here("output", "tables", "tab_loss_risk.txt")
modelsummary(
  m.risk,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefmap,
  gof_omit = "R2|IC|Log|F|RMSE|Err",
  add_rows = data.frame(
    c("Covariates"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes")
  ),
  gof_map = gof_map,
  fmt = 2,
  output = "latex"
) %>%
  add_header_above(c(" " = 1, "Donation" = 2, "Program" = 2, "Information" = 2)) %>%
  cat(., file = risk.file)
fix_txt(risk.file)


# Heterogeneous treatment effects by trust
m.trust <- list()
# Donation allocation outcome
m.trust[["Own"]] <- lm(update(f.base, scale(donateloss) ~ loss * index_trust + .), g.tot)
# Neighbors allocation outcome
m.trust[["Neighbor"]] <- lm(update(f.base, scale(neighbors) ~ loss * index_trust + .), g.tot)
# Interest in participating
m.trust[["Participate"]] <- lm(update(f.base, interest_num ~ loss * index_trust + .), g.tot)
# Support for program
m.trust[["Support"]] <- lm(update(f.base, losssupport_num ~ loss * index_trust + .), g.tot)
# Interest in receiving information
m.trust[["Seek"]] <- lm(update(f.base, info_num ~ loss * index_trust + .), g.tot)
# Share information
m.trust[["Share"]] <- lm(update(f.base, share_num ~ loss * index_trust + .), g.tot)

trust.file <- here("output", "tables", "tab_loss_trust.txt")
modelsummary(
  m.trust,
  vcov = "HC2",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = coefmap,
  gof_omit = "R2|IC|Log|F|RMSE|Err",
  add_rows = data.frame(
    c("Covariates"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes")
  ),
  gof_map = gof_map,
  fmt = 2,
  output = "latex"
) %>%
  add_header_above(c(" " = 1, "Donation" = 2, "Program" = 2, "Information" = 2)) %>%
  cat(., file = trust.file)
fix_txt(trust.file)
